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Abstract 

We investigate the use of stochastic methods for zero energy quantum scatter- 
ing based on a path integral approach. With the appHcation to the scattering 
of a projectile from a nuclear many body target in mind, we use the potential 
scattering of a particle as a test for the accuracy and efficiency of several meth- 
ods. To be able to deal with complex potentials, we introduce a path sampling 
action and a modified scattering observable. The approaches considered are 
the random walk, where the points of a path are sequentially generated, and 
the Langevin algorithm, which updates an entire path. Several improvements 
are investigated. A cluster algorithm for dealing with scattering problems is 
finally proposed, which shows the best accuracy and stability. 



1 Introduction 



The most frequently used method to study small scale composite systems is quantum 
scattering. While the interaction of the projectile with an isolated constituent of the 
target may be well known and can be described accurately, the scattering from the 
bound system can in general only be solved with several simplifying assumptions. 
These approximations concern both the interaction of the projectile with the target 
as well as the internal target dynamics during the scattering process. 

While a better understanding of the scattering from composite systems is of 
interest for several fields of physics, we focus here specifically on the scattering of 
a strongly interacting particle from nuclei at zero energy. There exist several often 
used recipes for the description of low enegy nuclear scattering, mainly in the form of 
effective single particle potentials for the projectile - nucleus interaction. However, 
the inherent approximations have not been tested in a conclusive way. That there are 
surprises if one goes beyond the standard approximations has been shown recently 
in the framework of the frozen nucleus picture, where all nucleons are held fixed 
during the scattering process. By evaluating the scattering amplitude with Monte 
Carlo methods for a target with a large number of scattering centers, it was found 
that the projectile gets 'trapped' by clusters of target constituents 0. This is an 
effect which is absent in any of the approximate descriptions. Clearly it is important 
to investigate methods how to solve in practice the scattering from a many-body 
target exactly, such that also the target dynamics can be taken into account. 

A formalism that has hardly ever been applied to scattering is the path inte- 
gral description. In a recent study ^ this approach was applied to the potential 
scattering of a particle and stochastic methods were used to evaluate the multi- 
dimensional integrals one encounters. It was then shown that these methods can in 
principle be extended in a straightforward fashion to calculate the exact scattering 
amplitude for scattering from nuclei with the full inclusion of the target dynamics. 
Before this method can be applied to a detailed study of nuclear scattering, which 
requires a considerable computational effort, it is necessary to find more efficient 
and numerically stable stochastic methods to obtain the scattering amplitude. This 
is the purpose of this paper where we examine different stochastic algorithms that 
can be apphed to the scattering problem. To be able to assess the accuracy of these 
methods we apply them to the potential scattering of a low energy particle, where 
the exact solution can be obtained by standard numerical methods from the single 
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particle Schrodinger equation. Keeping in mind the application to problems such as 
the interactions of low energy anti-protons with nuclei, we include complex poten- 
tials in our investigation of stochastic methods for scattering, which has not been 
considered before. 

In several fields of physics path integrals are a well established approach to ob- 
tain the ground state properties of a system and there exist stochastic methods 
to evaluate the resulting high- dimensional integrals. Examples are spin systems 
IQ, and quantum fields on the lattice 0. However, the subject of these stud- 
ies are localized, bound systems. The application of these methods to scattering, 
where the states are not normalizable and the spectrum is continous, is therefore 
not straightforward. Another complication is due to the use of complex potentials, 
which have been already used for bound states e.g. in quantum chemistry 0, |^ or 
in the calculation of correlation functions of nuclear ground states 0. A complex 
action makes it impossible to interpret the distribution of paths as a probability 
distribution and we examine how to deal with this problem in scattering. For com- 
pleteness we mention here also the possibility of obtaining the scattering length of 



a system from the volume dependence of its energy spectrum |10|. The feasibility 



of this approach has been shown for some special cases of quantum field theory on 



the lattice: the non-linear cr model in two dimensions |Tl| and the Ising model in 
four dimensions ||12|. The generalization of this method to nuclear scattering is cer- 
tainly not straightforward and we therefore focus below on more direct path integral 
methods. 

In Section 2, we briefly outline the path integral approach to non-relativistic 
quantum scattering and define a 'modified observable' appropriate for scattering 
calculations with a complex action. Different stochastic methods to evaluate the 
resulting multi-dimensional integrals are discussed in the subsequent chapters. We 
examine their efficiency and accuracy in order to identify the method that is most 
suitable for the application to many-body targets. The random walk algorithm, the 
most straightforward implementation of the path integral concept, is discussed in 
Section 3. Langevin and hybrid methods, which update entire paths, are studied in 
Section 4. An efficient new approach, a cluster algorithm for quantum mechanical 
scattering, is developed in Section 5. A summary and our conclusions are presented 
in Section 6. 
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2 Potential Scattering 



2.1 Scattering formalism 

In this chapter we briefly describe the theoretical framework for solving a scatter- 
ing problem by using the path integral approach. We discuss only the potential 
scattering for a single particle. It is used in this paper as a test case for finding a 
numerical method that is suitable for applications to scattering from a composite 
many-particle target. Some modifications of standard methods in the literature are 
necessary since we allow for complex potentials. 

The Hamiltonian, which will be used throughout this paper, describes the non- 
relativistic scattering of a particle from an complex, local potential 



H = ^ + V = k + V, 
2/i 



V 



U + iW, 



(1) 
(2) 



where p is the momentum operator and fi the reduced mass of projectile - target 
system. We are interested in stationary scattering states, corresponding to a real 
energy E{k) = Due to the complex potential the Hamiltonian H is not hermitian. 
Therefore one has two solutions to the Lippmann-Schwinger equations, which belong 
to H and W, respectively |T^: 
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We assume the existence of a complete biorthogonal basis for the Hamiltonian, H 

The scattering wave functions | for zero projectile energy can be obtained 
from the propagator /)(/?), 



J\l/3 L 



/27r/i^ 



(4) 
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in the limit of large Euclidean times (3. This operator contains the complete physical 
information about the evolution of the system. Its spectral representation in the 
biorthogonal basis given by Eq. is 



m = d'k exp [-PE{k)] 



In the limit f3 ^ oo one is left with 



lim p{P) 



(6) 



(7) 



which projects onto the ground state of the system. This projector can be used to 
calculate the square of the scattering amphtude, 

2 



(27r)V (Xo|^|^o'")(^o iV^IXo). 



(8) 



The wavefunctions | xo) ^-^e plane wave solutions, which in the zero energy limit 
are constant. The wave functions | tp^) and | ijj^) differ with respect to their 
boundary conditions. However, since no spatial direction is distinguished in the 
limit E{k) — > 0, Eq. (§) therefore can be written as 

2 



or in r-space representation 



(27r)V \im{xo\V p{P)V\xo), 



(9) 



(27r) /i lim / dridrN+i Xo{rN+i)V{rN+i)p{rN+i, (3) K(fi)xo(n) • 



(10) 

Note that, in contrast to other discussions of scattering in the literature, the quantity 
/q is calculated instead of | /o p. It contains information not only about the absolute 
value, but also the phase of the zero energy scattering amplitude. 

In the following chapters, we will investigate different stochastic methods to 
solve Eq. (p!oD . In all subsequent applications, the shape of V will be assumed to 
be Gaussian, 



V{f) = {Uo +tWo)v{r) = V,v{r) , 



v[r) 



exp 



262 



(11) 
(12) 



but this assumption is not crucial for the methods discussed here. Keeping in mind 
applications to the scattering length of antiprotons |T^ , |[T^ , we take the real part of 
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the potential Uq in the range of —0.5 fm ^ to 2.5 fm ^ and Wq from to — 7.0 fm ^. 
The width parameter will always be taken to be 6 = 0.5 fm and the reduced mass 
fi = 2.5 fm-^ 

2.2 Path integral approach and modified observable 

For the actual evaluation of Eq. we use the Feynman path integral approach. 
For this purpose the projector, Eq. (|^), is rewritten as [jT8| 



A[R] 



N 

P-l^e „=2 



N+1 

E 

n=2 



exp 



2vri 



/3 = eN. 



(13) 
(14) 

(15) 



In the integrand of the path integral A[R] denotes the Euclidean action and R = 
[fi,f2, ...,?>/] a path consisting of steps. Each step corresponds to an increment 
e of the imaginary time, /3. Parts of the potentials appearing in Eq. (|^) will be 
included into the projector by extending the action to 

A[R] = A[R] - Invi - InvN+i , fn = v{fn) , 



or in our case to 



A[R] = A[R] + 



262 



^2 I ^2 



Combining different factors into 



AT 



(16) 
(17) 

(18) 
(19) 

V[R] = dfidf2...dfN+i , (20) 

the expression for the scattering amplitude, Eq. (^, finally becomes a multi- 
dimensional integral 

f^ = ^JV[R]P4R]. (21) 



and with the definitions 



and 



27r/i2 ' 
Pa[R] = exp \-A[R] 
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The integral now also includes the endpoints fi and tat+i of the path. 

The expression for the scattering amplitude, Eq. (|2T|), contains only known 
quantities, e.g. it does not involve the unknown scattering wave function. Instead 
one has to deal with multi-dimensional path integrals. In the next chapters we 
will examine different stochastic methods to evaluate these integrals with sufficient 
accuracy and speed, such that the path integral approach can be also be applied to 
many-body targets. A problem is that the potential appearing in the action A is 
complex. The expression in Eq. (^1]) therefore cannot be interpreted as an integral 
over paths sampled according to a real and positive distribution, the starting point 
for many stochastic methods. Even for real potentials, which have both attractive 
and repulsive regions, this interpretation requires special care [§. 

To remedy this difficulty associated with complex potentials, we introduce a 
suitable real action, S, that yields a real positive distribution. 



Ps[R] 



exp 



-S[R] , Ps[R]>0. 



(22) 



The action S should simulate the main features of the actual problem as closely as 
possible or, equivalently, P^f-R] provide a good approximation to the kernel P4[i?] 
in Eq. (0). For a many-body target, for example, such a reference problem could 
be described in terms of a simple effective one body potential [Q: in regions where 
the actual interaction is strongly absorptive and the wavefunction is small, a real 
reference potential Ua{r) may be choosen that is strongly repulsive. Using this 
approximate reference action, Eq. ( pi] ) can be rewritten as 

'Pa[R] 



' c 



V[R] Pj,[R] 



c 



V[R] 



Ps[R] 



Ps[R] 



(23) 



The expression has the form of a multi-dimensional integral over a generalized com- 
plex observable, 0[-R], where the distribution of paths is given by Ps[R]'- 



f 2 _ 



0[R] 



V[R] 0[R] Ps[R] = 

Pa[R] 
Ps[R] 



{O) 



Ps 



or more explicitely 



jf / nR] e 



(A[R]-S[R])^-S[R] 



c 



^5 



(24) 
(25) 

(26) 
(27) 
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where 



AS[R] = A[R] - S[R] . (28) 

As we will see, in working with complex modified observables one encounters strong 
numerical fiuctuations which are caused by the presence of a phase that varies along 
the entire path. 

Monte Carlo methods used to generate paths distributed according to Ps always 
yield a normalized distribution. Except for few cases, this normalization is not 
known. Therefore, instead of Eq. (|26|) one calculates the quotient 



Ul !V[R]e-m UliOa)' 



(30) 



where Ua is the strength of the real reference potential. The advantage of this 
method is obvious from Eq. (|30D : one does not calculate the whole problem from 
the beginning, but uses a solved one as a reference problem. The computational 
effort is spent on calculating just the difference between both. 



8 



3 Random- Walk- Algorithms 



The most direct implementation of the path integral expression for the scattering 
amplitude is the method of random walks. We will first discuss this approach for real 
potentials and later apply it to complex potentials. The procedure for generating 
paths starts with a set of initial points ri. From each of these points the next point 
is then reached by means of the Euclidean propagator (r, | e""^ | ri). The set of 
paths generated in this way is then distributed according to Ps- This procedure is 
continued until a sufficiently large time f3 is reached to allow the projector e~^^ to 
filter out the zero energy contribution. As the projection time grows, more paths 
outside the potential region contribute, while configurations inside, which sample 



the dynamics of the target, only yield a small contribution [|T9|. This causes large 
fiuctuations in the results and we examine methods to improve the generation of 
paths through the guided random walk (GRW) algorithms. 

3.1 Simple methods: free diffusion and diffusion with sources 
and sinks 

Free diffusion: The simplest method to generate paths uses only the free action Sq 
in the Euclidean propagator, 

So[R] = i: ^{r^^i-rnf. (31) 

n=l 



The matrix element in Eq. (p3|) is then evaluated by using paths distributed accord- 
ing to 

Ps[R] = jj^^ exp [-S,[R] + Int^i] , (32) 
Mo = J dfv{f) < oo . (33) 



The real potential U only enters through the generalized observable 0[-R], 

K 



0[R] = 2^1? Ul ^ exp \\yivn+i - Su[R]] , (34) 



N+l 

Su[R] = eJ2 Un. (35) 

n=l 

The procedure to generate paths according to the distribution Ps[R] consists of the 
following steps: 
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1. Choose starting coordinates r/"\ a — 1, according to the initial 
distribution $^(ri), 

^.(n) = ^t'(ri) . (36) 

We only discuss here the case that v{f) does not change sign, which is true for 
our choice of a Gaussian shape. 

2. Each of the Zj„ points is successively propagated with the transition proba- 
blility p, 



p(4,ni+i;e) = -^exp 



9 . V n '^n+i; 



(37) 



for the step from point n to n + 1 of a given path a. This is numerically 
implemented by taking 



ri?! = ri"^ + y'^ r/i"^ , a = l,2,..Z,„, (38) 

where is a vector consisting of three independent, uniformly distributed 
random numbers. 

3. Step 2 is repeated N — ^ times. 

4. Finally, the scattering amplitude is obtained as 

■^tn a=l 

Replication of paths: Properties of the potential can be incorporated into the 
generation of paths by the method of sources and sinks for paths: each of the 
segments of a path is assigned a weight Dl^\ 

Di-) = exp[-£C/(r<'^))]. (40) 

After each free diffusion step as described above, the segments of the path generated 
up to this 'time' ne are replicated or deleted, depending on the value of D^""^: 

Let Dint be the largest integer which satisfies < Dint < and draw a 
random number w in the interval [0,1]. For w < (D^°^) - Ant) one makes 
Dint copies of the path segment up to this point, otherwise only {Dint — 1) 
copies. These additional paths are then continued independently. If however 
Dint — 0, the segment will be kept with probability w or deleted with the 
complementary probability {1 — w). 
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It is easily seen that this leads to the replication of paths in the region of attractive 
potentials; for repulsive potentials there is no replication of paths and the methods 
ensures that paths with low weight factor D^"^ are not followed too far. Especially for 
the case of strongly attractive potentials, not all paths one obtains in this fashion are 
uncorellated, since the paths from replication contain common segments. Compared 
to the free diffusion, the advantage of this method is that the paths are distributed 
according to Ps, where S is the full action, i.e. they are selected in a way that 
reflects the dynamics of the system. 

We now illustrate the above by applying these simple random walk methods to 
a repulsive and an attractive real potential. To be able to judge the quality of the 
methods, we calculated the exact expectation value of the observable as a function 



of the Euclidean time (3. This was done by solving Eq. (|T^) in its phase space 
representation with the fast Fourier transform (FFT) method [^. Figs. |l| and |^ 
show results obtained with the free diffusion and the path replication method. By 
comparing to the FFT calculation, we see that for an attractive potential the con- 
vergence is much slower and relatively large projection times, typically f3 = 100/m, 
are needed. This long projection time requires high stability of the calculations. In 
the following we will therefore often use the attractive potential as a test case. 

The difference in quality between the two methods, free diffusion and path repli- 
cation, is small. The main source of this insensitivity can be traced to the rapid 
diffusion from the potential region,which leads to a relatively small number of paths 
that actually contributes to the matrix element of the observable: starting with 
100000 paths, the number grows to 135000 for the attractive potential. Of these, 
only 630 contribute at /5 = 100/m, i.e. the generation of new paths does only little 
to balance the loss of paths due to the free diffusion from the potential region. For 
a repulsive potential, the number of paths drops from the initial 100000 to 50000, 
of which only 120 contribute. For free diffusion, the total number of paths stays 
constant and 330 contribute. The number of contributing paths is not sufficiently 
different in the cases considered to lead to a noticable difference in accuracy. In view 
of the achieved accuracy and the long projection time, none of the methods seems 
suitable for application to complex obsevables and many-body targets, which have 
higher demands on numerical accuracy and speed of the calculations. 



11 



3.2 Guided random walk algorithms 

An improvement of the random walk method can be achieved if an approximate 
solution for the wavefunction of the system is known. This solution can be used for 
the generation of the paths one uses to sample a modified observable. For the scat- 
tering from a many-body target such an approximate solution can be obtained for 
example from an effective single particle potential Of course, since we consider 
the potential scattering of a particle, the exact wavefunction in our present stud- 
ies can be obtained by solving the Schrodinger equation with standard numerical 
methods. To get an impression of the power of the guided random walk method, we 
use here the 'optimal' situation: we work with the exact wavefunction for the real 
potential. This provides an upper limit for the performance of this method, since 
with only an approximate wavefunction the accuracy can be expected to get worse. 

For the zero energy case, we thus use the solution of the Hamiltonian with a real 
potential, U{r) = Uov{f), that satisfies 

H ^(f) = 0. (41) 

The 'time-independent guided random walk method' only makes use of this station- 
ary solution. We have tested the time independent method to generate paths and 
evaluate the appropriate observable. It was found that in our case this method also 
could not prevent the substantial drifting of paths out of the potential region. Con- 
sequently, there was no noticeable improvement in accuracy compared to the simple 
random walk method. We therefore describe here only the time dependent guided 
random walk, which tries to compensate the loss of configurations. An essentially 
constant density of configurations in the potential region can be ensured by using a 
zero energy wavefunction with a norm that depends on the Euclidean time, Te, 

One then considers the evolution of the product (j){'r,Tf.)ip{'f,Tf.), which has the in- 
finitesimal evolution operator 

Ps = <p{f,T, + e) e~'^ <P'\f,Te) . (43) 

It can be shown PT| that the matrix element can be written in terms of integrals 
over paths distributed according to a modified distribution 

' ^ ( 1 [ n ^ ^ ^^ 

v{ro) , (44) 



2 En /i 
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V In ^„ , €n = e 



1 + -A ln< 



and a simple modified observable, 



0[R] 



/i(2 7r/i) 2 Uq Afy exp[\nvN+i - In^jv+i + In^i 



After discretization, the infinitesimal evolution operator reads 



Pe 



n + 1 



n 



AC. 



exp 



2 e. 



[fn+l - r„ Fnf 



(45) 



(46) 



(47) 



The numerical implementation of the above is similiar to the diffusion case: 
successive points of a path are obtained according to 



' n+1 



n II 



Vn 



(a) 



where a = 1, 2, ..Z,-, 



After each step, there is a replication of the path segment by a factor 



n + 1 



n 



(48) 



(49) 



through a procedure analogous to the one in Section |3.1| . The difference is that now 
the weight factor is not space-, but time-dependent, reflecting the time-dependent 
norm of 0. This change of the particle number was found to be successful in keeping 
the number of points in the potential region at a reasonably high level. The result 
is an increase in accuracy and stability, as shown in Fig. ^ 

The price one pays is that the total number of points grows rapidly with time. To 
keep about 2000 points in the potential region over a projection time of /3 = 100/m, 
an initial number of 100 configurations grows to ca. 200000. Since the generated 
paths all go back to a very small number of initial configurations, there is a noticeable 
correlation among the paths for small (3, which can be seen by the large initial 
deviations from the exact result. However, for large projection times the method is 
remarkably stable. This is also due to the presence of the drift-force F in Eq. (^), 
which together with the random configurations t] has a decorrelating effect. Such 
effects will be discussed in more detail in the next chapter in conection with the 
Langevin algorithm. 

Due to the reasonable stability, the time-dependent GRW method looks like 
a promising method for the complex potentials, which require higher numerical 
accuracy. The full Hamilton-operator H is now 



H = K + U + iW . 



(50) 



13 



The wavefunction 0(r, t^) is determined through the real potential with the time- 
independent part satisfying 



[K + U{f)] e(r) = . (51) 

The modified observable also includes also the absorptive part of the potential 
through the imaginary action Sw [R] , 

0[R] = fi{2nfi)-^U^ K exp[lnf^+i - In^^+i - Sw[R] + In^i]. (52) 

The resulting scattering length is shown in Fig. ^. The time dependent GRW does 
somewhat better than the random walk method based on the free action. The figure 
shows that there are less fiuctuations in the results obtained with the time-dependent 
GRW. However, some fiuctuations remain. They are mainly due to the phase, which 
sums up the contributions along each path. Similar results were also found for an 
absorptive potential with a repulsive real part and are not shown. 

In summary, we have seen that the random walk algorithm converged reasonably 
fast towards the exact scattering length. However, the numerical accuracy was not 
satisfactory for any of the different versions we tried. A typical calcuation took ca. 
25 minutes on a HP UX 9000/720. Thus while these algorithms are not useful for 
the calculation of scattering from a many-body target, they can provide estimates 
for the typical projection times one needs. These projection times are needed as an 
input in the algorithms described in the next chapters, where paths are not built up 
successively from point to point, but entire paths of a given length /? are generated 
and modified. 
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4 Langevin and Hybrid Algorithms 



4.1 Standard Langevin Algorithm 

The path integral methods presented in the last section describe a diffusion process 
in the physical time /3, sequentially generating the elements of a path, Another 
approach is the Langevin algorithm, which updates an entire given path R. This 
development of a path proceeds in terms of a new variable r, the 'Langevin time', 
and takes place according to the stochastic differential Langevin equation, 

|-i?,(r) =-S,'iT) + n,{T) , (53) 

where vr is a Gaussian random variable 

(7r,(r)) =0 , (7r,(r),7r,(r')) = 2S,,S{t-t') . (54) 

The variable r is different from the physical time, which is represented here by the 
index i. When discretized, r labels the elements in a Markov chain. Again, iS is a 
real approximation to the action A and 

SAr) = ^mr)]. (55) 
The stochastic process described by Eq. ( |53D is ergodic, 

{0)p ^ f V[R] Ps[R] 0[R] = Jim i r dr' O [R{t% (56) 



where 



Ps\R] = ^e^^\-S[R]] . (57) 



s 



Choosing for S the real part of the action, the paths generated through Eq. ( p3D 
will be used to evaluate the modified observable 

0[R] = exp [-5vk[^]1 . (58) 



For numerical evaluation, Eq. (|5^) is discretized. We first write the contributions 
to the action S, 

S[R]=So[R]+Su[R], (59) 
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in a convenient matrix notation 

So[R] = So[R] + ^ (fi^ + rt+i) = RAoR, (60) 

N+l 

Su[R] = eJ2U, = -TrDu[i?] , (61) 
i=i ^ 

^ f lUv{f^ forz = l,iV + l 

\ t/t;(fi) for i = 2,3,...,iV . ^ ^ 

The matrix Aq is symmetric and has 3 (A^ + 1) eigenvalues, Aj > 0, while Du[-R] is 
diagonal and built from N blocks, 

Bu[R] = [Ui Es] , 1 = 1, 2, {N + l) , (63) 

where E3 denotes the three-dimensional unit matrix. With these definitions we 
obtain the discretized version of the Langevin equation that updates a path i?*^") 
and yields R^''+'^\ 

^("+^) = [1-6 M(")] + h , (64) 
M(") = Ao-^Du^"^ . (65) 

The continous variable r in the Langevin equation has been replaced by the discrete 
label n, 

n = ^ , (66) 
where 5 is the stepsize and we have defined h through 

(67) 

The Langevin algorithm thus generates from a starting path a sufficiently 
large number of path updates through Eq. (0), which are used to evaluate the 
observable. However, not all successive paths produced by the iteration of Eq. (|5^) 
are independent and to assess the efficiency of the algorithm we have to consider the 
autocorrelation time T^or separating two decorrelated paths. As can be seen from 
Eq. (0), the evolution in r is governed by the eigenvalues /ij of the matrix M. It 



can be shown that 

Tcor ~ , (6^ 
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where /imm is the smallest eigenvalue of M. On the other hand, for the algorithm 
to be numerically stable, the step size 6 should be small enough to resolve also the 
evolution of the quickly changing modes, i.e. 

6 « . (69) 

Consequently, the number of intermediate steps ricor, required such that a path i?*^™) 
is decorrelated from its predecessor R^'^\ satisfies 

Tcor fJ'max . /^n\ 

ricoT = -7- >> , m = n + n^or- (70) 

In general, the spectrum of the full matrix M cannot be calculated. Since our 
scattering potential has a short range, the potential matrix Du^"'' vanishes in most 
regions of the space of paths. As an approximation, we therefore replace the full 
dynamical matrix, M, by the matrix containing the free action, Aq, and get a lower 
limit for the autocorrelation time in terms of the eigenvalues Aj, 

^ ^ ^max \ 
ricor » T • (71) 

The ratio of the largest and smallest eigenvalues is shown in Fig. 5. As can be 
seen, a large number of intermediate steps must be taken before an uncorrelated 
path is obtained. This is due to the presence of two very different length scales in 
scattering: The average the distance of the middle of the path from the origin is 
of order 0[jjj and grows with the projection time. On the other hand, since the 
endpoints of the paths are required to lie in the potential, a small scale b enters. Of 
course, the actual spectrum of eigenvalues in scattering is continous and the presence 
of a lowest eigenvalue A^m > is due to discretization. The finite projection time, 
the existence of a condition for the endpoints and the discretization of the action are 
thus crucial for the convergence of the Langevin algorithm in scattering problems 

i- 

The change of the stability criteria due to the potential can be roughly estimated 
in the region where Du is approximately constant. For a path entirely in this region 
one has 

^ « ^^2) 

(,2 

For an attractive potential, numerical stability may therefore require a stepsize 5 
which is considerably smaller than the estimate based on the eigenvalues of Aq; 
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for repulsive potentials the stepsize may be taken larger, at least as long as the 
values of Uq are such that the r.h.s. of Eq. ([7^) stays positive and finite and the 
condition doesn't become meaningless. The presence of the potential thus may 
require an increase in the number of intermediate steps one needs before another 
uncorrelated path is obtained or, even worse, it can lead to to numerical instabilities. 
The condition Eq. ( [7^ ) does not apply in the region where the potential varies 
significantly. Lowering the stepsize in the physical time, e, to reduce the influence of 
the potential on the stability it is no remedy, since the eigenvalues Aj also dependend 
on e. However, since the range of the potential is relatively small, most steps take 
place in the outer region and Eq. (|68| ) still provides a reasonable indicator for the 
overall autocorrelation time. 

Finally, to avoid local instabilities and errors due to the discretization of r, 
the updated path obtained from the discretized Langevin equations are treated as 



the proposal step in a Metropolis procedure p6|. The resulting algorithm is often 
referred as the 'exact' Langevin algorithm. 

The above considerations are illustrated by the examples in Figs. |] and ^ For 
a potential with a strong repulsive real part and strong absorption as in Fig. |^, a 
projection time of /? = 50 fm and e = 0.5 fm is sufficient; as in Section this was 
checked by comparing to the results of deterministic FFT methods [^]. Inspection 
of Fig. 1^ shows that one needs according to Eq. (|7lD Ucor » 5000 intermediate 
Langevin steps in r to generate an independent path. The numerical values for 
Xmax and Xmin are 20 and 4.7 ■ 10^^. The stepsize 6 is chosen to be 2.5 ■ 10~^ in 
order to meet condition Eq. (^) at least in the region where the potential vanishes. 
Therefore the result shown in Fig. ^ contains less then 80 decorrelated paths. After 
an initial equilibration phase of about 5000 Langevin steps, the the fluctuations get 
smaller and the result for the real part of the scattering length gets very close to the 
exact value. However, the imaginary part is still quite far off and shows little sign 
of further convergence to the exact answer with increasing r. 

For attractive potentials the situation is worse and the necessary projection time 
increases. For the strength parameters in Fig. |^, Uq = — 0.3fm~^ and Wq = 
— l.Ofm"^, the FFT calculations show that a projection time P = 80 fm should 
be chosen. This introduces a larger length scale and the number of intermediate 
steps Ucor increases quickly. For e = 0.5 fm the number ricor is already of the 
order of 10000 and less than 40 independent paths contribute to the result in Fig. 
^ Correspondingly, one has Xmax = 20 and A^m = 1-9 ■ 10"'^; 6 is chosen to be 



18 



1.0 • 10~^. The equilibration time is now larger, ca. 10000 steps, and fluctuations in 
the result remain as r increases. One improvement would be to reduce discretization 
errors in (3 for the longer projection time by reducing the stepsize e. This would lead 
to an even larg illustrated by Fig. ^. However, the calculation in Fig. |^ 

already took approximately 2.5 hours on one processor of a CONVEX EXEMPLAR 
1000 using standard tools for algebraic manipulations. Therefore this possibility of 
improving the performance seems not practical. 

The rather poor performance of the Langevin algorithm seen in the two examples 
considered here is very different from the successful application to real potentials in 
Ref. P|, which could even be extended to many-body targets. The major cause of 
this difference is the phase in the modifled observable, Eq. (^8|), which is responsible 
for strong fluctuations. Clearly the few independent paths in the examples above 
are not sufficient to probe the rapidly varying observable and their number has to be 
increased to achieve a satisfactory result. Straightforward extensions of the above 
method by increasing the number of Langevin steps would lead to prohibitively long 
computation times. We therefore now consider two possibilities which can help to 
accelerate the algorithm: The first one modifies the dynamical properties of the 
Langevin process by introducing a kernel into the stochastic differential equation, 
Eq. (|53D, while the second method, the hybrid algorithm, combines the Langevin 
algorithm with a microcanonical approach. 



4.2 Langevin Algorithm with Kernel 



To speed up the convergence of the Langevin algorithm, an appropriate matrix 
kernel, K, can be inserted into the discretized Langevin equation, Eq. (|64|). This 
does not change the distribution -Psf-R] according to which the paths are sampled. 
The conditions under which the use of a kernel is valid are discussed in e.g. Refs. 



P2| - p4|. For position- independent kernels the equation becomes 



2 



Ao 



+ /i Tic 7f . 



(73) 



Note that the square root of the matrix K re-scales /i, the stepsize in the 'Langvin 
time' T. The dynamical matrix, M, which now governs the evolution of the algorithm 
in r, is 



K 



Ao 



(74) 



with n 
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Again the region outside the potential dominates the dynamical properties of 
the algorithm. To increase the efficiency, it is desirable to decouple the slow, long 
ranged modes of the free action So, corresponding to the small eigenvalues of Aq, 
from the fast modes. A natural choice for the kernel is K = Ao~^, which leads to 

M(") = 1-^[Ao-^DuH, (75) 

reducing to the unity operator outside the potential region and thus has eigenvalues 
Aj = 1. This lowers the correlation time, Eq. (^), quite dramatically to Tcor ~ 1! 

As discussed in the previous section, numerical instabilities can occur in the 
region of the potential. With the choice K = Aq^^ the potential does not enter 
into the kernel and the scaling of the stepsize in the inner and the outer region of 
the potential is the same. Nevertheless, since the correlation time is dominated by 
the behavior in the outer region, where the different eigenmodes are decoupled by 
the use of the above kernel, far less intermediate steps are required to generate the 
necessary independent paths. 

The improvement of the algorithm through the kernel does not permit an increase 
in the stepsize 6. The coupling of the modes in the interior remains and the numerical 
stability requires again to follow each step proposed by the Langevin equation with 
a kernel by a Metropolis acceptance-step. Since the algorithm is dominated by the 
long ranged modes of the free action, a large rejection rate is obtained. This has to 
be remedied by a decrease in the step-size 6 compared to the case without kernel 
in order to have the Metropolis step accept 50% to 80% of the proposed paths. 
Another factor that leads to longer computation times is the matrix structure of the 
kernel and thus the increased algebraic complexity of the algorithm. The number of 
the required numerical operations scales at least quadratically with the dimension 
of the matrix Aq. 

Fig. P shows the result for a complex potential with a repulsive real part. The 
calculation with the projection time (3 = 80 fm, consumed approximately 2 hours, 
the runs with /? = 50 fm and different discretization e took 1 hour and 3 hours, re- 
spectively. The computational effort and the accuracy is similar as in the calculation 
without a kernel in Fig. ^ But the increased number of about 5000 independent 
paths leads to a very stable behaviour; furthermore the number of steps for the 
initial equilibration is now reduced to about 1000. As can be concluded from Fig. ^, 
a smaller discretization in the projection time (3 is needed to improve the accuracy 
of the result. This trend was veryfied by a run with e = 0.1 fm, /? = 50 fm and 
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n = 5000 Langevin steps, which however took more than 5 hours (not shown). 

Fig. ^ shows results for a potential with an attractive real part and /3 = 80 fm, 
which took about 2 hours. Although the results are again remarkable stable after 
relatively few Langevin steps due to the use of a kernel, they deviate from the exact 
results. The calculations were repeated (see Fig. [1^) with a longer projection time. 



(3 = 200 fm, and fewer Langevin steps, n = 5000, which required 5 hours of CPU 
time. The results are somewhat improved, but are less stable. The fluctuations 
are caused by the decreased number of independent paths. Calculations for more 
attractive potentials were also carried out and it was found that no reasonable 
convergence could be achieved within the chosen time limits. 

The usefulness of introducing a kernel into the Langevin algorithm is therefore 
limited to a small range of potential parameters, where it can produce very stable re- 
sults. However, as this improved algorithm is rather time-consuming, its application 
to multi-particle systems looks not promising. 

4.3 Hybrid- Algorithm 

The Langevin algorithm can be combined with deterministic molecular-dynamical 
methods to a 'hybrid algorithm' with the advantages of both methods 
PU| . The deterministic updating of entire paths in a new 'hybrid time', t, is governed 
by a 'CPU Hamiltonian', Ti, which is different from the physical Hamiltonian, H, 

n = \n' + S[R] (76) 

- (77) 

dt^' - ^^^^ 



d on 



dt OR. 



-S'm- (78) 



In this approach vf is the CPU momentum of an object with position vector i?, mov- 
ing in a potential S[R\. The discretized version of these equation for the transition 
from n to n + 1 reads to order 0[/i^] 

After several of deterministic steps according to the above equations of motion, the 
momentum vf is updated by a Gaussian-distributed random number. Therefore, all 
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momenta and position vectors, if they they belong to a deterministic step or not, 
are random variables. 

— * 

The trajectory of a path R is split into sub-trajectories, 7^, each starting with a 
Gaussian random momentum, jf('^i°'\ and ending with vf^"-^'*), containing 

Aria = — nia- (80) 

sucessive deterministic steps T2 (/;,). The next sub-trajectory 7^+i then starts again 
with a Gaussian random momentum. Each move along a sub-trajectory is com- 
pleted by a Metropolis acceptance step: The entire sub-trajectory is considered as 
a proposed move, which is accepted with probability 



Parr. = miu 



^11 



AH = 7^[^("^"\7f("^")] -7^[^("^"\7f(''^'^)] . (82) 

This combination with the Metropolis step, which we have used in all examples 
below, is referred to as 'exact hybrid algorithm'. 

The length of 7^, i.e. the number of deterministic steps Ana, should be such that 
the paths and are decorrelated. Note that the first two steps of Eqs. 

(|79D can be recombined to yield the ordinary Langevin equation, Eq. (0). In the 
Langevin case, a sub-trajectory T then just consists of one step and the Metropolis 
step in both methods is analogous. The Hamiltonian Ti. describes a classical motion 
confined in a high- dimensional harmonic oscillator potential. In contrast to the 
previous approaches, the correlation between two paths both belonging to the same 
sub-trajectory has not an exponential but oscillatory dependence on the time which 
labels successive paths. The period of this dependence is ti = where Aj again 
denotes the eigenvalues of the free action matrix Aq. To assure decorrelation of the 
smallest mode within one sub-trajectory, its length should satisfy 

1 

Ura ~ ry^ ■ (83) 

However, this does not prevent strong correlations between some other modes due to 
their different periodic time dependence of their correlation. Therefore in practice 
the length of each sub-trajectory is varied randomly in the vicinity of ttra- Stability 
requires a number, ritra, of steps in a subtrajectory that satisfies 



tfra ^max /r>A\ 

= -7- ~ \/t • (84) 
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Comparison of Eqs. ( ^l)) and (^) shows that the decorrelation in the hybrid algo- 



rithm is faster by a factor i/if^ compared to the standard Langevin algorithms; the 
substantial size of this factor can be seen in Fig. ^. In addition to the decrease in the 
required number of intermediate steps, a further advantage of the hybrid algorithm 
is that it needs fewer Gaussian random numbers, which are very time consuming to 
generate. 

Eqs. (|79D represent only the simplest discretization in the time t. With increas- 
ing sub-trajectory length the discretization error increases also the drift from the 
CPU energy shell in the deterministic steps of a sub-trajectory. This drift is mea- 
sured by AH, Eq. (^21) . Therefore the acceptance rate is lowered considerably and 
an improvement of the discretization in Eqs. (^) is necessary. We have studied 
examples with projection times in the range /3 = 50 - 200 fm and stepsizes e = 0.1 
- 1.0 fm. For sub-trajectories containing approximately ritra = 50 steps, the Olh"^] 
discretization was found adequate. For sub-trajectories of about 100 steps a dis- 
cretization procedure up to terms of order 0[h^] and for 200 steps even 0[h^] was 
needed to lower A7i, requiring three or nine times more steps, respectively, as with 
the 0[/i^] discretization. We used the higher order scheme of Ref. [^, where Eq. 



(|7|) is extended to 0[h''+^] by the recursive relation 



Tk+2{h) = Tk{h)Tk{-sh)Tk{h) (85) 
h = T^^h, s = 2^. (86) 

We found that in the range of parameters shown in Fig. ^ the hybrid algorithm, 
while technically more involved, works efficiently; discretization to order 0[h^] was 
sufficient. 

The calculations shown in Figs. |ll] and are obtained with the exact hybrid 
algorithm, without a kernel. The potential is Vc = (2.3, —7.0) fm~^, the projection 
time (3 and the discretization e are varied. As before, we found for this repulsive 
potential a projection time of /? = 50 fm sufficient. In this case discretization in the 
stepsize, h, to order 0[h^] was needed for e = 0.33 fm, while for £ = 0.1 fm also terms 
of order 0[h^] were necessary. The calculations of discretization, e = 0.33 fm, and 
order, 0[/i'^], consumed about 30 minutes of computing time, while the e = 0.10 fm, 
0[h^] example needed more than 5 hours! Note that compared to 0[h^] calculations 
not only the number of algebraic manipulations is increased by a factor of three, 
but according to Eq. ( p^ and Fig. |^, longer projection times (3 and smaller values 



23 



of e also lead to an significant increase in ritra- The length of one trajectory can 
be estimated to range from ttra ~ 15 to ttra ~ 30 (see Fig. H). This means, that 
Fig. contains approximately 6000 independent samples for the discretization 
e = 0.33 fm and 2500 samples for e = 0.10 fm. The algorithm shows a good and 
stable performance. Reasonable results can be obtained for a smaller number of 
steps in the 'hybrid time' t which required about 10 minutes of computational time. 
Figs. |lT] and ^ show that the discretization is important, but a slight discrepancy 
in the imaginary part still persists. 

The limitations of the method show up for attractive potentials. While the 
hybrid algorithm was found to work equally well for potentials with weak attractive 
parts down to —0.3 fm~^, it starts to fail for more attractive values of Uq- This is due 
to the fact that with our standard parametrization of 6 = 0.5 fm and yU = 2.5 fm~^, 
a real potential develops a bound state at Uq = — 1.2fm^^. Fig. |13| shows the the 
results for Vc = (— 0.5, — 1.5) fm~^. The time needed for the calculations was 90 
minutes for a projection time f3 = 100 fm. Attempts with (3 = 200 fm (not shown) 
were also found insufficient to filter the ground state well enough and therefore the 
algorithm cannot be used for attractive potentials of this strength. 

Both extensions of the Langevin algorithm — use of a kernel and hybrid method 
— show a comparable behaviour. In our present applications of the algorithm, all 
versions comprise a sequence of steps given by Eq. ([79|); if a kernel used, h is scaled 
by the square root of K. For stability reasons h was chosen to be approximately five 
times smaller then minimum required. The hybrid algorithm performs slightly faster 
because it needs far less Gaussian random numbers. For refined discretizations and 
long projection times, the two improved versions become equally time-consuming. 
It is hard to compare both extended schemes in general, since the calculation time 
depends on the algebraic structure of the chosen kernel, the implementation of the 
stochastic and deterministic differential equations and the discretization in h for 
the hybrid algorithm, which will vary from case to case. A kernel can also be 
introduced into the hybrid algorithm. We did not examine this possibility in detail, 
but tests showed that long trajectories cause stability problems. Therefore, one has 
to decrease the step size further and no gain in velocity is obtained. 
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5 Cluster- Algorithms 



The hybrid and Langevin Algorithms are hindered by the fact that the couphng 
of two points of a path is restricted to nearest neighbor interactions. The physical 
potential is purely local and acts only within a small region with typical scale b of the 
potential width. Outside this region one is left with the free action. Hence, strongly 
oscillating paths are suppressed and long range correlations of order /? dominate. 
Only the end point contributions containing the logarithm of the potential prevent 
the paths from drifting away to infinity. Hence, the displacement of a path on a 
large scale proceeds very slowly. 

To improve the computational efficiency it is desirable to speed up these slow 
long-ranged modes without neglecting small fluctuations introduced by the short 
range of the potential. One would like a non-local or collective-mode updating 
based on information over a larger range of coordinates. We considered already one 
example of a such an algorithm in Section ^ where we introduced a kernel into the 
Langevin algorithm in order to deal with modes of the path motion and not just 
individual points when generating new paths. This approach suffered from the fact 
that we were not able to diagonalize the complete action. A coupling of the different 
modes through the potential remained. 

A similiar situation often occurs for spin systems with nearest neighbor interac- 
tions. It is well known that large-scale modes of spin waves lead to sizeable slowing 
down. To remedy this, Swendsen and Young introduced a new type of collective 
algorithm for Potts models . Additional auxiliary degrees of freedom are intro- 
duced into the original model: groups of correlated spins are treated as new degrees 
of freedom which are called 'clusters'. Changes are applied to these collective objects 
and not to single physical degrees of freedom, the indivual spins, which leads to a 
much more efficient and faster way to deal with such systems. Subsequently many 
related applications have been developed for Ising models, 0(N) spin models p and 
lattice gauge theories I^Tj, commonly referred to as 'cluster algorithms'. There is 
considerable freedom in how to define a cluster. Some special aspect of the physical 
problem has to be used in each case to define the appropriate collective variables. 

In this section, we extend the idea of clusters from bound systems to scattering 
problems. Adjacent coordinate points of a path having similiar positions will under 
certain conditions be considered as clusters and dealt with collectively. The size 
of these clusters must refiect the two different length scales of the physical system 
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under consideration, the short range of the potential and the large exterior region. 
There are two basic steps of the cluster algorithm which are carried out to a given 
path R to generate a new path, B! . They are repeated until satisfactory convergence 
of the matrix element of the chosen obserevable is achieved: 

1. The first step consists of grouping the points of a path into clusters. As will 
be discussed in detail below, clusters are designated in a way that will favor 
a decrease in kinetic energy when updating the path by moving these clusters 
in step 2. 

2. Once the clusters have been defined, a number of them are reflected on a 
randomly oriented plane yielding a new path E! . This new path consists 
thus of some clusters which were reflected and some which did not change. 
The decision wether a cluster is reflected or not is chosen in way way which 
favors the decrease of the potential energy. To preserve ergodicity, the position 
of £ is changed each time this step is carried out. 

Both steps together lead to a splitting into smaller clusters in regions where the 
scale of the potential becomes important and to larger clusters where the influence 
of the kinetic part of the action is dominant. This is a great advantage compared 
to e.g. the methods in Section ^ There are two differences between our method 
and the application of cluster methods to spins of fixed length. First, it is in our 
case not sufficient to reflect the coordinate points through a plane which contains 
the origin; instead the position of the plane is varied throuhgout the entire space. 
Second, the scattering potential influences the reflection probability of clusters and 
therefore scales the effective size of reflected clusters. 

5.1 Introduction of auxiliary variables and collective updat- 
ing of paths 

To group points of a given path R into clusters, first auxiliary 'bond' variables hi are 
defined for each pair of adjacent points and rj_i. The hi can take on the values 
or 1. The whole set of these variables for a path is denoted by 5. A cluster C is a 
set of neigbouring coordinates fi which are connected by non-zero bonds hi 



C = [ri,ri+i 
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where s stands for the set of indices i,i + 1, j of coordinate points contained in 

To decide which value should be assigned to a 'bond' variable 6j, we use a plane 
S, characterized by a normal vector n and its distance ^ to the origin. We then 
consider the point f*, which is obtained by reflecting fi through the plane S, i.e. 
f* = f*[fi,S]. The probablitiy of setting the bond variable to zero is then given 
by 



^nb,i = mm 



l,exp(-A5o[A,]) ] , (88) 



where ASq is determined by the free part of the action, 

A5o[A,] = |i[Ar-A,2], (89) 

A, = f^-f^-l, (90) 

A^ = r7-rl_i. (91) 
P[B I R, S] is then conditional probability to obtain a set of bonds B for given plane 

N+l 

P[B \R,8] = n \Pnh,i Sb,,o + (1 - Pnb,) , (92) 

1=2 

B = [b2,h, ....,bN+i] R= [ri,r2,----,rN+i] ■ (93) 

Taking into account these additional variables, the probability distribution Ps [R] 
can be extended to the combined distribution V 

V[R,S,B] = P[B I R,S]Ps[R] Pe[n,^] . (94) 

Here Ps[R] is as defined before in Eq. (|57D and Pslh,^] describes the probabihty of 
finding a plane S with normal vector n and at a distance to the origin. Its form 
will be specified below; in general it does not depend on the path R. One has by 
construction 

Ps[R] = I dij d^Y. V[R,£,B] , (95) 

{6} 

where integration over d Q extends over all possible orientations of n. In terms of the 
extended probality distribution, V[R,S, B], we can re-write the expectation value 
of the modified observable 0[R] as 

(O) = I V[R] 0[R] Ps[R] = f V[R] j di j dVtY ^[R] V[R,S,B] . (96) 
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We now have to specify how we update the path R by reflecting entire clusters 
on the plane S. We denote by the cluster one obtains by reflecting C'^ on The 
probability for a cluster to be reflected is given by 

T[C'*\C'] = amin[l,exp(-A5t/[C"]) ] , (97) 
ASu[C'] = Su[C'']-Su[C']. (98) 

Here < a < 1 is a real parameter; its choice will be discussed below. The quantity 
5[/[C'^] is the contribution of the potential to the action for the points within the 
cluster C*. Clusters which contain the endpoints also contribute Inf according to 
Eq. ([16|) . As can be seen from the above prescription, the probablity for reflecting 
a cluster is chosen such that a transition to clusters with a lower potential energy is 
favored. Since the new path R' in general contains reflected clusters, C**, as well as 
unreflected ones, C*, the total transition probability for given S and B reads 

TsbIR' \R] = l[ T[C'* I n {1 - ^P"* (99) 

s=l n=l 

The first factor in the above equation contains only the clusters which are re- 
flected to transform R into R', while the second is due to the A^^^ clusters which 
retain their position. After having updated a path R to R', the procedure is then 
repeated with now R' as starting path to obtain a new R' etc., until a sufficient 
number of paths is obtained for the evaluation of the observable. 

It remains to be shown that the transition from R to R' satisfies the criterion of 
detailed balance, 

nR',£,B] ^ Tsb\R' I R] .^QQ. 

V[R,£,B] TsbIR \R']' 
Making use of the general relation 

min[l,e~''] _^ 



min[l, e' 

one finds for the r.h.s. of Eq. ( |100| ) 

Tsb[R' I R] 



(101) 



exp 



s=l 



(102) 



TsbIR I R'] 

The sum runs only over those A^,. clusters which are reflected to make the transition 
from R to R'. The l.h.s. of Eq. ( |100| ) can be re-written by making use of Eq. ( |101| ) 
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to yield 



V[R,S,B] 
where the exponent is given by 



exp 



-AS[R', R] 



P[B I R',£] 
P[B \R,£y 



(103) 



Ntr Nr 

AS[R',R] = S[R']-S[R]=J2ASo[A,] + J2ASu[C' 

t=l s=l 



(104) 



The index t in Eq. ([ToiD runs only over those Ntr junctions which are changed 
into A*; junctions which are contained in one cluster and junctions between two 
clusters which are reflected simultaneously do not contribute. The remaining term 
in Eq. ( p.03| ) is obtained from Eqs. (|88D , (p2D, which yield 



P[B I R,£] 
P[B I R',£] 



exp 



Ntr 

t=i 



(105) 



Combining Eqs. ( |103| ), ( |104| ) and ( |105| ) and inserting the result into the l.h.s. in 
Eq. ( |100[ ), it is readily seen that the detailed balance condition is satisfied. 



5.2 Scaling of the Cluster Size and Choice of Parameters 

Important for the convergence and accuracy of the algorithm are the size and relative 
distribution of the clusters. Too many large clusters leads to inefficient probing of 
the potential region, while too many small clusters slow down the algorithm. The 
cluster distribution depends on the function P^ and the coefficient a in Eq. (p7|). 
We will discuss these two aspects in this section. 

Geometrical considerations show that two neighbouring points belong to different 
clusters if the line connecting these points is intersected by the plane £. Furthermore, 
it is readily seen that the size of the path segment contained in a cluster and the 
relative distribution of clusters with different sizes depend on the average distance 
of the plane £ to the region of classical paths where there are many points of a path. 
A distribution P^-, which mainly places the plane into these physically preferred 
areas, results in many intersections and thus a large number of small clusters. If 
on the other hand the plane £ is placed at large distances from these regions, it is 
likely to intersect the path only at a few points, resulting in relatively few, large 
clusters. Thus, while at first sight the setting of bonds and grouping of points into 
clusters seems to entirely depend on the kinetic energy through the condition in Eq. 
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(|88D , the above considerations show that also the potential enters into the eventual 
distribution of clusters through the action S: If the algorithm is equilibrated and 
starts to converge, the paths will be denser in regions of local minima of the real 
action S; reflection planes placed into these regions will lead to a splitting into many 
small clusters that allow for an accurate probing of the details of the potential. 

The location of the plane S is important for another reason. If S is placed be- 
tween two classically allowed regions, which are separated by a repulsive potential 
well, the reflection makes the fast interchange of path segments between both min- 
ima possible. A very large number of intermediate steps would be needed if, instead 
of clusters, we dealt with individual points as in the other algorithms. With the 
Langevin method, for example, the path flrst has to get through the repulsive poten- 
tial barrier step by step before it reaches the neighboring minimum. This advantage 
of the cluster method is thus especially advantageous for scattering problems with 
a more complicated target structure such as many-body targets. Projectile paths 
will not remain trapped in regions separated by potential barriers, but move quickly 
and probe efficiently the entire target region. 

We now discuss the choice of the parameters defining the distribution P^ln,^] 
of the refiection planes. The orientation of the plane S is defined by the normal 
vector n; in spherical coordinates n = {cosip sin ^, sin?/' sin 0, cos 0). Since we are 
considering scattering problems for spherically symmetric potentials no direction 
in space should be preferred and therefore the variables cos 6 and ip are randomly 
distributed over the intervals 

cos^Gi-l,!], ij e[0,27r]. (106) 

The distance of the plane to the origin, ^, is taken to be of Gaussian type. 



1 



1 

2 

exp 



2w 



' e 



(107) 



If the projection time is large enough such that it satisfies << -, it can be shown 



that one has for the mean square distance of the midpoint of the paths to the origin, 

(4.) ^ 3|^, (108) 
while the mean square distance of the endpoints is 

«^+i)^3 6^ (109) 
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For repulsive potentials the choice w ~ 6 for the width of the distribution leads to 
the best results. By concentrating the position of the reflecting plane in this region 
around the origin, one avoids that large segments in the exterior region, typically 
located at a distance {(Pm) » b are reflected into the potential region. Avoiding 
these paths, which are practically all rejected, speeds up the algorithm and reduces 
fluctuations. At the same time, this positioning of the reflecting planes takes into 
account the detailed shape of the potential in updating a path. On the other hand, 
taking w too small would lead to too many almost symmetrical reflections, as the 
planes are concentrated close to the symmetry center of the potential. The updated 
paths then only differ in the kinetic energy, while the contribution from the potential 
remains essentially unchanged. As the observable we evaluate only depends on the 
absolute value of the distance of the points of a path to the origin, such updates 
would only increase the computational time for the algorithm. In the examples 
below with a repulsive real part of the potential, which were all carried out with a 
projection time oi P — 80 fm, we chose a width oi w — 0.5 fm. 

Attractive potentials require a much longer projection time. As a result, the 
midpoint of the paths moves further away from the origin and a larger distribution 
parameter w ~ 0[dM/2] should be chosen, yielding planes £ located about halfway 
between the midpoint of a path and the origin. Exterior segments are then often 
reflected into the energetically favored potential region; the average time a cluster 
remains in the interior region is short, since the next path- update will tend to move 
it out of that region again. Symmetric reflections through the origin are suppressed 
with this choice for A disadvantage of this choice results for the endpoints of the 
path, which should remain in the potential region: a large w leads to reflections out 
of this region and thus to conflgurations with a low acceptance. We used e.g. for a 
projection time oi P — 300 fm a width oi w — 2.0 fm, while for P — 500 fm a choice 
oi w ^ 4.0 fm yielded the best results. 

The above discussion shows that one has to balance competing effects when 
choosing the distribution of the planes: If w is chosen too small the tendency exists 
that many clusters are reflected to a position at the same distance to the origin; a 
too large w leads to motions of the endpoints with a low acceptance. Both extremes 
signiflcantly lower the convergence speed of the algorithm. Clearly a different shape 
of the surface on which one reflects is also possible to avoid these problems and to 
further improve the efficiency of the algorithm; we do not pursue this here. 

The probabihty for the reflection of a cluster depends on the potential and on 
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the choice of the parameter a in Eq. (0), which also has an important influence on 
the relative distribution of clusters of different sizes. The decision if a cluster should 
be reflected is made independently for each cluster. If we take a ~ 1, Eq. (p^) 
shows that a single cluster is then practically always reflected if the potential energy 
is not increased by this operation. Therefore in regions with no structure due to the 
potential, it is then likely that many adjacent clusters are reflected simultaneously 
and 'effective clusters' of large size are formed. Within the range of the potential 
small clusters are formed since the reflection probability, Eq. (pTj), now changes 
over short distances due to changes of the potential. The algorithm then provides 
a fast collective movement in regions where the potential vanishes and leads at the 
same time to a fine scaling when the path crosses the relatively small interaction 
region. Therefore the choice of a ~ 1 is especially suited for scattering problems. On 
the other hand for a = 1 a cluster would always be reflected and effective clusters 
would tend to stay together in the exterior. This leads to a slowing down due to 
insufficient decorrelation of modes with a longer range. We obtained the best results 
with a 0.75. In choosing the parameters to optimize the cluster motions in the 
exterior we do not have to worry about how this affects the modified observable. In 
our case the observable depends on the imaginary part of the potential and has no 
structure outside of the range of the potential. The choice of a can thus entirely be 
used to speed up the algorithm in the large exterior region. 

Fig. |l^ demonstrates for a repulsive potential how the algorithm selects auto- 
matically the appropriate scale for the clusters. Small clusters, necessary to probe 
the potential region in detail, can be found at short distances from the origin, while 
larger clusters can be found at larger distances. That the cluster size does not grow 
stronger at very large distances from the origin is due to the fact that the endpoints 
of the path are restricted to lie in the potential region. Too large clusters lead to 
too few acceptable paths. 



5.3 Examples 



We first show in Fig. 15 results for a potential with a repulsive real part. The 



parameters are the same as in Fig. |TT], where it was shown that with the hybrid 
method a discretization of the projection times /3 down to e = 0.25 fm could not 
achieve convergence to the correct answer in both real and imaginary part. The 
results were obtained with the cluster algorithm for a refined discretization of e = 
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0.1 fm. There is now a rapid convergence to the correct result for the scattering 
amphtude. For the hybrid algorithm such a fine discretization would have been 
extremely time consuming and out of the question for an extension to many-body 
targets. In contrast, the calculations in Fig. ^ took less than 60 minutes with the 
cluster method. (All calculations discussed in this section were preceeded by 50000 
equilibration steps before starting with the evaluation of the expectation value of 
the modified observable; the quoted computation times always include this initial 
equilibration.) That the refined discretization was essential for the accurate result 
in Fig. |15| can be seen in Fig. |l^, where different choices for e are compared. 



Comparison of Figs. |Tl] and |I6| shows that the influence of the discretization of the 
projection time is equal for both methods; a calculation with the same projection 
time (3 and discretization e yields results of comparable quality for both the hybrid 



method used in Fig. |TT| and the cluster method. 

To right away also test the limits of the cluster method, we used attractive 
potentials which, as we have seen before, put higher demands on the stochastic 



algorithms. The same qualitative results as for the repulsive potential in Figs. 15 



and were also found for weakly attractive potentials as long as the strength does 
not exceed Uq = —0.3 fm~^. However, for potentials with a stronger attractive part 



the necessary projection time increases rapidly. This can be seen in Fig. |T^, which 
shows the dependence of the accuracy on (3. With the efficient cluster method one 
can afford to go to projection times as high as 900/m to achieve good agreement 
with the exact result; a coarser discretization than for the repulsive potential could 
be used. The computation time with the cluster method was about 120 minutes for 
the (3 = 900 fm case. Calculations with such a long projection time are impossible 
with the methods of Section § within the limits we set for computation time. The 
results shown in Fig. |l3|for the same attractive potential and with a projection time 
(3 = 200 fm required 180 minutes and show no signs of convergence to the exact 
result. 



Closer inspection of Fig. [T^ shows that even for the best result some fluctuation 



persist. This feature is even more evident in Fig. |1^, where the attractive part of 
—0.75 fm~^ required a very long projection time, resulting in a computation time 
of 180 minutes. As discussed above, for high projection times the middle of the 
path moves out of the potential region. Reflections occur increasingly on planes at 
large distances from the origin, causing the endpoints of the path to move out of 
the potential region. This leads a low acceptance, slowing down the speed at which 
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uncorrelated paths are generated and fluctuations remain in the result. 

To give a general impression of the performance of the cluster algorithm, we 
show in Fig. [TP| a set of results for attractive and repulsive potentials with varying 
imaginary parts. The corresponding parameters are listed in Table |l]. After the 
equilibration, 300 000 new paths were generated in all cases to evaluate the observ- 
able. The error-bars in Fig. |T9| are obtained from the results of the last 50 000 
steps. The solid curve shows the exact result obtained by numerically solving the 
Schrodinger equation with a Runge-Kutta method. The computation time for the 
repulsive potentials was 45 minutes, while the attractive potentials require a longer 
projection time f3 and took ca. 60 minutes each. The overall performance is quite 
good. The only significant discrepancies occur for the attractive potential, which is 
close to having a bound state and requires longer projection times. The examples 
with a weak imaginary potential could actually have been carried out in about half 
the computing time, since the scattering wave function of the full problem is quite 
similar to that with the real part of the potential only (which is used to generate 
the paths) and the convergence is fast. 

The strength of the real simulation potentials we have used here varied from 
Uq = —0.75 fm~^ to 2.3 fm~^. To put this range into perspective, we show in Fig. 
^ the corresponding exact scattering length as a function of the strength parameter 
Uq. The scattering length changes substantially and grows rapidly as the real part 
of the scattering potential approaches Uq = —1.2 fm~^ where a bound state appears 
and /o diverges. 

In summary, the performance of the cluster algorithm is superior to that of any of 
the other algorithms discussed before. While the hybrid algorithms can be improved 
by finer discretization, the necessary computation time grows roughly quadratically 
with the number of points |jl9| and becomes unreasonably large. In contrast, the 
increase in computation time is only linear for the cluster algorithm, which allows one 
to obtain results quite efficiently even for the difficult cases of attractive potentials 
which the hybrid algorithm cannot manage anymore. In cases where both hybrid 
and cluster methods could be applied, the achieved accuracy was comparable. To 
further improve the succesful cluster method and understand its performance in 
more detail, one would need to better understand the decorrelation of paths in this 
method. Further studies should also be directed towards optimizing the orientation 
of the reflection plane to yield a uniform acceptance of reflections of both the middle 
segments and end points of a path. 
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6 Summary 



The path integral approach offers an opportunity to solve exactly the quantum me- 
chanical scattering of a hadron from a composite, many-body target. The standard 
methods, which have been used for decades, involve several approximations con- 
cerning both the reaction mechanism between projectile and target as well as the 
internal target dynamics. Recent studies have shown that due these approxima- 
tions one misses important physics, such as the trapping of the projectile by several 
scattering centers. 

The possibility to apply the path integral formalism to low energy scattering 
has been pointed out several years ago. However, little has been done since then to 
actually use it in practice and to investigate which computational method is best 
suited for its application. The path integral approach is of course a common tool for 
the study of localized, bound systems. But when considering scattering one is faced 
with new aspects. An obvious difference is that in scattering the wavefunction is not 
normalizable and the spectrum of the Hamiltonian continous. Another important 
feature that has to be dealt with by the numerical methods is the presence of two 
very different length scales. One is provided by the size of the target or the relatively 
short range of the interaction between projectile and the target constituents. The 
other scale, which is much larger, is introduced by the exterior region into which 
the paths can drift off with increasing physical time. 

In order to assess the efficiency of the different methods for applying the path 
integral approach to scattering, we have tested them for a situation where the exact 
answer can still be obtained by standard methods: the scattering of a particle from 
a potential. As the idea is to find a method that can be extended to the much 
more extensive and time consuming calculations for a many-body target, the desired 
computational method has to be fast and numerically stable. Keeping in mind the 
applications to systems with absorptive channels, e.g. for calculating the antiproton 
- nucleus scattering length, we have used a complex potential. Due to the presence 
of a complex action, the usual probability interpretation of the distribution of paths 
is not possible. We therefore used the fact that one can choose a different real 
action to generate the paths and to accordingly rearrange the matrix-element of the 
observable one evaluates. In our approaches, the imaginary part of the potential is 
always contained in such a 'modified observable'. In the action used for generating 
paths we used here for simplicity the real part of the scattering potential; in other 
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applications it is of course possible to make another choice suggested by the physics 
of the system. 

We compared three different methods to generate paths for the evaluation of the 
observable. The first is the random walk method, where paths are built up point 
for point. In contrast, the Langevin algorithm, with its extensions to the Langevin 
algorithm with a kernel and the hybrid algorithm, deals with the updating of entire 
paths, rather than stepping through each path in the physical time. These two 
methods, which had previously been applied to bound systems, showed systematic 
difficulties when apphed to scattering. We could trace the difficulties of these meth- 
ods, which have been used extensively for bound states, to the different physical 
situation in scattering and tested some modifications and hybrid methods. The in- 
sights gained from these studies finally this led us to propose the 'cluster algorithm' 
for scattering. It deals not with individual points of a path but instead with clusters 
of varying sizes. 

For the random walk method we considered both path generation by the simple 
free diffusion method as well as the method of path rephcation. Both yielded rather 
inaccurate results, mainly due to the fact that only few of the generated paths 
actually contributed significantly to the integral. Some improvement was obtained 
by the 'time-dependent guided random walk', but large fluctuations in the observable 
persisted, in particular for attractive potentials, which require a large projection time 
and put higher demands on numerical stability. 

While the random walk sequentially generates the elements of a path, the Langevin 
algorithm updates an entire path. As discretization errors become important, the 
generated paths were treated as proposals for a Metropohs procedure. Within 
the computational limits we considered, this procedure was seen to be quite time- 
consuming and yielded only unsatisfactory results. This was traced in part to the 
small number of uncorrelated paths produced in this fashion. To speed up the algo- 
rithm, a kernel was introduced into the stochastic differential equation, which allows 
one to decouple the long-ranged and short-ranged modes and substantially lowers 
the 'decorrelation time' between independent paths. This yielded very stable results, 
but only over a limited range of potential parameters. As the matrix algebra makes 
this method more time consuming, its application to more complicated scattering 
problems does not look promising. A similar performance was found when combin- 
ing the Langevin method with deterministic methods to the 'hybrid algorithm'. This 
was shown to lead to a substantial gain in the the decorrelation time. Nevertheless, 
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the need for refined discretization also makes this method rather time consuming, 
especially when attractive potentials are considered. 

Instead of dealing with individual points when updating a path, the 'cluster 
algorithm', which we finally proposed, treats groups of points as the new degrees of 
freedom. When generating a new path from a given one, these clusters are reflected 
collectively on a randomly oriented plane. Due to the procedure for forming and 
reflecting these clusters, this was seen to lead to an updating in terms of large scale 
clusters in the exterior region, while small scale clusters are used in the potential 
region, allowing an accurate probing of the details of the potential. This algorithm 
was found to be very stable and efficient, able to handle all but the very strongly 
attractive potentials. Another aspect that makes this method look very promising 
for application to nuclear targets is the reflection procedure. It allows a path to move 
quickly over larger distances, e.g. between regions that are separated by repulsive 
potential barriers. If point-by-point updating were used, this would be extremely 
time consuming. 

For none of the three types of algorithms have we tried to optimize specific 
computational aspects, such as the matrix manipulations for the Langevin algorithm 
with a kernel. Individual improvements are clearly possible and will cut down the 
computation time in each case. However, we believe that the qualitative differences 
among the different approaches which we have found are quite general and that the 
cluster approach is the one that should be pursued first when considering scattering 
from nuclear targets. 
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Tables and Figures 
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Table 1: Potential parameters (in fm ^) used for Fig. 19. 
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Figure 1: Scattering length /o for repulsive potential, Uq = 
1.0 fm~^, as function of projection time, /3; discretization e — 
0.5 fm. Open squares: with path replication; solid squares: 
free action random walk; solid curve: exact FFT result. 
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Figure 2: Scattering length /o for attractive potential, Uq — 
—0.3 fm~^, as function of projection time, /3; discretization 
e — 0.5 fm. Labelling as in Fig. 1 
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Figure 3: Scattering length /o for attractive potential, 
Uo — —0.3 fm~^, as function of projection time, /3; discretiza- 
tion e — 0.5 fm. Open squares: time dependent GRW; solid 
squares: free action random walk; solid curve: exact FFT 
result. 
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Figure 4: Scattering length /o for complex potential Vc — 
(—0.3, —1.0) fm~^ as function of projection time P; discretiza- 
tion e = 0.5 fm. Open squares: time dependent GRW results; 
solid squares: free action random walk; solid curve: exact 
FFT result. 
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Figure 5: Upper graph: smallest eigenvalues of free action ma- 
trix Aq as function of discretization e for different projection 
times (5. Lower graph: ratio of largest and smallest eigenvalue 

of Aq. 





Figure 6: Scattering length for complex potential Vc — 
(2.3, —7.0) fm~^, obtained from the 'exact' Langevin algo- 
rithm; (3 — 50 fm, e — 0.33 fm. Horizontal line: exact result. 





Figure 7: Scattering length for complex potential Vc — 
(—0.3, —1.0) fm~^, obtained from the 'exact' Langevin algo- 
rithm; P — 80 fm, e — 0.5 fm. Horizontal line: exact result. 
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Figure 8: Scattering length for complex potential Vc = 
(2.3, —7.0) fm~^, obtained from the 'exact' Langevin algo- 
rithm with kernel. Horizontal line: exact result. 
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Figure 9: Scattering length for complex potential Vc — 
(—0.3, —1.0) fm~^, obtained from the 'exact' Langevin algo- 
rithm with kernel; P — 80 fm, e — 0.5 fm. Horizontal line: 
exact result. 
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Figure 10: Scattering length for complex potential Vc — 
(—0.3, —1.0) fm~^, obtained from the 'exact' Langevin algo- 
rithm with kernel; P — 200 fm, e — 0.5 fm. Horizontal line: 
exact result. 
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Figure 11: Scattering length for complex potential Vc = 
(2.3, —7.0) fm~^, obtained from the 'exact' hybrid algorithm. 
Horizontal line: exact result. 
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Figure 12: Scattering length for complex potential Vc = 
(2.3, —7.0) fm~^, obtained from the 'exact' hybrid algorithm. 
Horizontal line: exact result. 





Figure 13: Scattering length for complex potential Vc — 
(—0.5, —1.5) fm~^, obtained from the 'exact' hybrid algo- 
rithm; P — 100 fm, e — 0.33 fm. Horizontal line: exact 
result. 
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Figure 14: Mean square distance to the origin as function 
of effective cluster size for /3 = 80fm, £ = O.lOfm, Vc — 
(2.3, -7.0) fm"^ . 



-0.850 



-0.900 




-0.950 

0.0 100000.0 200000.0 300000.0 



0.40 




E 



0.20 
0.15 

0.10 

0.0 100000.0 200000.0 300000.0 



Figure 15: Scattering length /o as function of the number 
of paths, Np, generated by the cluster algorithm for Vc — 
(2.3,-7.0) fm~\ /3 = 80 fm, £ = 0.10 fm. Solid line: exact 
result. 
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Figure 16: Scattering length /o as function of the number 
of paths, Np, generated by the cluster algorithm for Vc — 
(2.3, —7.0) fm~^, P — 80 im and different discretizations, e. 
Solid hne: exact result. 





Figure 17: Scattering length /o as function of the number 
of paths, Np, generated by the cluster algorithm for Vc — 
(—0.5,-1.5) fm~^, different discretizations, £, and projection 
times, p. Solid line: exact result. 
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Figure 18: Scattering length /o as function of the number 
of paths, Np, generated by the cluster algorithm for Vc — 
(-0.75,-2.25) im-\ P = 900 fm, e = 0.50 fm. Solid line: 
exact result. 
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Figure 19: Scattering length for complex potentials with dif- 
ferent absorptive parts given in Table 1; absorptive part grows 
in absolute value from right to left. Upper figure: attractive 
real part; lower figure: repulsive real part. 
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Figure 20: Scattering length /o as function of strength pa- 
rameter Uq. Range of real simulation potentials used for path 
generation indicated by diamonds. 



